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A new algorithm is given for computing the solution of any second-order linear difference equation 
which is applicable when simple recurrence procedures cannot be used because of instability. Com- 
pared with the well-known Miller algorithm tin- now method has the advantages of (i) automatically 
determining the correct number of recurrence steps, (ii) applying to inhomogeneous difference equa- 
tions, (iii) enabling more powerful error analyses to be constructed. 

The method is illustrated by numerical computations, including error analyses, oi Anger- Weber, 
Struve, and Bessel functions, and the solution of a differential equation in Chebyshev scries. 
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1. Introduction 

A powerful computational algorithm for evaluating the most rapidly decreasing solution of a 
second-order homogeneous linear difference equation was published in 1952 by J. C. P. Miller 
(fl|J page xvii) in connection with the tabulation of modified Bessel functions. Since then, various 
writers have applied the algorithm to other special functions, and similar computational processes 
have been used by Clenshaw [2] for the numerical solution of ordinary differential equations in 
series of Chebyshev polynomials. Error analyses of the algorithm have been supplied by the present 
writer [3] and Oliver [12] and quite recently Gautschi |4J has examined the relation of the algorithm 
to classical results in the theory of continued fractions. 

The present investigation stems from the observation that Miller's algorithm can be regarded 
as a procedure for solving a tridiagonal set of simultaneous linear algebraic equations. Adopting 
this more general standpoint, we shall show how to recast the algorithm into a new form which 
enables the correct number of recurrence steps to be determined automatically without appeal to 
an asymptotic or other analytical formula. In this respect it resembles an algorithm proposed 
recently by Shintani [5]. 

The new formulation has the further advantages of (i) being applicable to inhomogeneous dif- 
ference equations, (ii) lending itself readily to powerful error analyses. There seems to be no 
alternative method of comparable power available at present for computing solutions of inhomo- 
geneous equations in the case when forward recurrence and backward recurrence are both un- 
stable. 

2. Statement of the Problem 

Let the given difference equation be denoted by 

a r y r -\ — h r Jr + C r y r + l = d r , (2.01) 



'Figures in brackets indicate the literature references at the end of this paper. 
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where a r , b r , c r , and d r are given functions of the nonnegative integer variable r. We assume that 
the general solution of (2.01) has the form 

yr = Af r + Bgr+hr, (2.02) 

in which A and B are arbitrary constants, and the complementary functions f r , g r , and the particular 
solution h r have the properties /o # 0, g, -^0 for all sufficiently large r, and 

frlgr^O, hr/gr^O, (r~*0o). (2.03) 

(It may be noted that we do not require either f r or h r to tend to zero as r^> oo.) 

The first problem we investigate is the computation of the solution of (2.01) which has the 
property 

Yrlgr^O (r-*oo), (2.04) 

and satisfies the normalizing condition 

yo = k (2.05) 

for an arbitrarily assigned value of the constant k. Later (sees. 9-11) we allow for a more general 
form of normalizing condition and also drop the restriction / ^ 0. 

The given conditions ensure that y r exists and is unique. For, from (2.03) and (2.04) the B of 
(2.02) is seen to be zero, and from (2.05) we derive A = (k — ho)lfo. Therefore 

k — ho r , , 

fr+hr. (2.06) 



Jt r 

h 

It is well known that direct use of (2.01) as a recurrence relation for generating y 2 , y.3, . . . 
from given values of yo and j\ (if available) is an unstable procedure. Essentially, each computa- 
tional rounding error introduces into the numerical solution a small multiple of f r and a small 
multiple of g r , and in consequence of (2.04) the latter ultimately grows faster than the wanted 
solution. 

It may also happen 2 in the inhomogeneous case that/ r grows more rapidly than y r in the direc- 
tion of decreasing r. In this event recurrence by use of (2.01) is unstable in this direction too. 

3. Approach 

Analogous work in the numerical solution of linear differential equations 3 suggests that a 
stable way of solving the present problem is to treat it directly as a boundary-value problem rather 
than use initial-value techniques. We are already given the value of y . Suppose that for some large 
integer TV, the value of y.v can be obtained from an asymptotic formula or by other means. Then 
eqs (2.01) with r= 1, 2, . . ., N— 1 comprise a set of simultaneous linear algebraic equations for 
the unknowns yi, y 2 , . . ., yv-i, which are solvable by standard matrix computational processes. 

This possibility has already been noted by Gautschi ([4], Introduction) following a suggestion 
by M. E. Rose, but Gautschi did not pursue the idea because of the difficulty of obtaining the value 
of y.v, in general. Following Miller's approach in the homogeneous case [1], we solve the algebraic 
equations with the value of y.v arbitrarily set equal to zero. It transpires that for large N the great 
majority of the y, produced in this way are generally excellent approximations to the true values; 
only in the neighborhood of r = N can substantial errors occur. 

We first establish the convergence of the process. 



2 See Example Lof section 6. 

3 See, for example, 16]. 
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THEOREM 1. With the conditions of section 2, suppose that for all sufficiently large N the sys- 
tem of equations 



aryffi - b r y r N » + c r y< r N + >, = d r , (r = 1 , 2 N - 1), 



(3.01) 



and 



has a solution y\^\ y*, 



N) v(N) 



y< N N) . 77^n i/r is fixed and N -> oo, y<N) -> y r . 



(3.02) 



To establish this result, we observe that since y r is a particular solution of (2.01), we can 
express 

W^Axfr + Bsgr + yr* 

where As and Bs are independent of r. Setting r = 0, A" in turn and using (2.05) and (3.02), we 
derive 



As = 



soys 



fogff-gofN* 



Bs 



/oy.v 



fogs-gofs" 



the denominator fog\~gof\ being nonvanishing for all sufficiently large /V in consequence of the 
assumed conditions. In fact, from (2.03) we have 



fogs-gofs~fo{ 



(N- 



>)■ 



Hence from (2.04) it follows that As ~* and Bs ~* 0. This completes the proof. 

Thus for any given value of r, or for any finite range of values, y r can be calculated to pre- 
scribed accuracy by solving the system of equations (3.01) and (3.02) with a sufficiently high value 
of N. Naturally, we now enquire what exactly constitutes a "sufficiently high" value? Or, to put 
the question another way, given /V, to what accuracy does ^ approximate y r for r < N? 

A simple practical way of providing an answer is to solve eqs (3.01) and (3.02) for increasing 
values of TV until the results are in satisfactory numerical agreement. This procedure has two 
drawbacks. First, it is wasteful of computing time if the originally guessed values of TV are either 
too low or much too high. Second, there is no absolute guarantee that values of y r computed with 
two (or more) different values of A^ must be correct when they agree. 

The optimum value of N, that is, the minimum value necessary to achieve specified accuracy 
in y r for a given range of values of r, can be determined automatically when a suitable method is 
used to solve the algebraic equations. Accordingly, we consider this process next. 

4. Solution of the Algebraic Equations 



We shall solve the tridiagonal system of equations 

a 2 ^-b 2 ^ + c 2 ^ = d 2 , 



a.v-2^3 - bs-2M 2 + CS-2M1 = d.V-2, 
a.v-i>^i-6.v-0^i=ifc-i, 



(4.01) 
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by simple elimination followed by back-substitution. To begin with we suppose that none of the 
c, vanish. 

Let the first of (4.01) be rewritten in the form 

p 2 ^)- Pxy ^) = e u (4.02) 



. b\ ci\k — d\ 

Pi = l, p 2 = ~, ei = (4.03) 

C\ Ci 

The result of eliminating y^ from (4.02) and the second of (4.01) can be expressed as 

where 

b-zp> — a 2 pi a 2 e\ — dkjk 

p.* = , ei = 

c 2 c 2 

Continuing the elimination, we obtain 

Pr+i&>- !*&$! = * (r=l,2,. . .,yV-l), (4.04) 

where 

b r p r — a r p r -\ a r e r - i — drPr 

Pr+i = ■> e r = (4.05) 



Thus p, is the solution of the homogeneous form of the difference eq (2.01), with the initial condi- 
tions p<> = and p\ = 1. We also observe that the second of (4.05) holds for r= 1 if we define 

e {) = k. 

The final equation of the form (4.04) is used to begin the back-substitution. On substituting 
the second of (3.02), we derive 

^-i^^v-i/pv; (4.06) 

thence y^l-,, ^-.^ • • ., y\ N) may be computed by use of (4.04) with descending values of r. The 
process fails if, and only if, one of the numbers p-2, P3, . . ., p,v vanishes. In this event the set of 
eqs (3.01) and (3.02) has either no solution or an infinity of solutions, and the algorithm breaks 
down. 

When one or more of the coefficients c r vanishes the set of eqs (4.01) becomes uncoupled. 
A simple modification takes care of the situation. Suppose, for example, that c s = but all other 
c r are nonzero. Then the first s equations of (4.01) determine y^, y{f\ . . ., y^ completely: they 
can be solved by application of the recurrence relations (4.05) for r= 1,2,. . .,5 — 1 and use of 
the back-substitution relation (4.04), beginning with 

y*)= SKt. (407) 

DsPs — dsPs-! 
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The remaining N— s— 1 equations are solvable for 3/^} , ^ , . . ., y^ by the method already 
described: eq (4.07) takes the place of the first of (3.02). 

To ease the presentation we shall suppose in the remainder of the paper that none of the c r 
vanish. 

Applying ourselves to the problem of determining the optimum value of /V, we observe that 
the effect of replacing /V by /V+ 1 is to prolong the elimination process by one step, beginning the 
back-substitution with yj^ l 1) =0 instead of y^ )= =0. Thus we have 

Pr + y r N+l) -pry^=e r (r= 1,2,. . ., TV). (4 .()8) 

Subtraction of (4.04) from (4.08) gives 



Pr+\ 



(4.09) 



and repeated application of this result leads to 



y, + „_y /; v, = J^££±i . . .£fizl (^V+D-y/)), 

Pr+lPr+2 PN 

that is, 

tf+»-flO = -&S!L- (r=l,2, . . ., N). (4.10) 

P.\P.\+\ 

By use of this formula we can predict the effect of changing N into /VH- 1 before any back-substi- 
tution is carried out. 

Suppose, for example, that we wish to compute yi to D decimal places for given values of the 
integers L and D. Then the recurrence relations (4.05) are applied from r=] past r=L until a 
value of r is reached for which 



Pi,e r 

PrPr+\ 



<JXHH>- (4.11) 



If this value of r is taken as N, then we can be sure that the approximation y ; v) yielded by the 
back-substitution agrees to D decimal places with the value y\ N+1) that would be obtained from 
the next higher approximation. (Whether this value of N is adequate is considered in the next 
section.) 

If, as is more usual, accurate values of y r are required for a whole range of values of r, then 
the criterion (4.11) is used with \pi] denoting the greatest value of \p r \ in the given range. We 
might, for example, desire the computation to I) decimal places of all values of y r that exceed 
V2 X 10 _/; in absolute value — it being assumed, of course, in this case that y r — > as r^> 00. Then 
/V is determined by the condition that 

kv/p.v+i|< JX10-", (4.12) 

provided that |pjv|^|pr| when r < N. 

5. Expansions for the Solution and the Truncation Error 

The method suggested in the last section for determining TV is based on the criterion that the 
values of y^ and )4 iV+1) must agree to within the prescribed tolerance for y r . This does not guar- 
antee, however, that their common value is y r . To resolve this doubt we consider higher approxi- 
mations. 
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Replacing TV by 7V+ 1 in (4.10), and adding the result to (4.10) itself, we obtain 

\PnPn+ i Pn+iPn+2/ 
Continuation of this process yields 



y* + ,)_y r iv) =Pr (_^+_^±I_ + . . . +. c ^-> 



PnPn+i Pn+iPn+2 Pn+s-iPn+s, 

where s is an arbitrary positive integer. Letting 5 — > °° and using Theorem 1, we derive the following 
expression for the truncation error 

gOv) = yr _ y n) = ENPr9 (5 01) 

where E\ is the sum of the (necessarily convergent) series 

E "=$ I T7- (5-02) 

Thus the precise criterion for determining TV is that \E\p r \ must not exceed the specified tolerance 
in y r for each wanted value of r. 

Once the value of N has been decided, the actual value of the truncation error can be found 
by continuing the computation of p r and e r -\ beyond r = N and using (5.01) and (5.02). Later [7], 
we shall show how to use these expansions to determine strict bounds for e^ directly from the 
properties of the coefficients a r , b r , c r , and d r . 

As a special case of (5.01) we have the expansion 

fT r PsPs+l 
Subtraction of (5.01) from (5.03) yields 

W^PrY-^- (r<N), (5.04) 

£ r PsPs+l 

a result which is obtainable more directly by repeated use of the back-substitution relation (4.04). 
Thus the whole of our computing scheme is equivalent to approximating the convergent infinite 
series (5.03) by the partial sum (5.04). 

6. Examples 

EXAMPLE 1. Anger-Weber functions. 

For integer values of r the function E r (^;) satisfies eq (2.01) with 

i a - 2r a - 2 n-(-i) r } 

a r =C r =l, O r — — , d r — 

X 7TX 

We restrict ourselves here to positive values of the argument x. The principal properties of E r (x) 
are established in [8], chapter 10. In particular, we have 

IM*) =— ^ | a,(r )x *, E 2r+1 (,) = <^ J ^^ a.(r) X >; 

(6.01) 
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where a (r) = 1, and 



a 8 (r) 



(4r*-&)(4r*-&) . . . Hr 2 - (2s + l) 2 } 



(s>0). 



(6.02) 



Using the inequality 



<**(r) 



we deduce that if x is fixed and r— > oo, then 

2* 



4r-l 



E 2 rU) 



(4^-1)77' 



E>2r+lU) ~7" 



(2r+l)77 



(6.03) 



The corresponding homogeneous form of (2.01) has the Bessel functions J r (x) and Y r (x) as 
solutions. For fixed x and large r, we have 



Mx) 



i 



(27rr) 1 / 2 \2r 



Y r (x) 



,_(2\ 1/2 (2r\ r 
\irr) \ex 



(6.04) 



Thus ultimately J r {x) decays more rapidly than E r (jc), and |Fr(#) | grows rapidly. In consequence, 
both simple forward recurrence and simple backward recurrence are unstable methods for gen- 
erating E r (x) from (2.01) when r > x. 
With 

fr = Jr(x), gr=Y r (x), kr=E r (x) , 

the conditions of section 2 are satisfied, provided that Jo(x) ¥" 0. Let us apply the method of section 
4 to a specific example, say the computation of E r (je) for x= 1, r= 1(1)10, correct to within 2 units 
of the eighth decimal place. We suppose that we are given E (l) = — 0.56865 6627, this value hav- 
ing been extracted from [9] and confirmed by evaluation of the first of (6.01). 

Beginning with p = 0, pi = l, and eo = Eo(l), values of p r and e r were generated by use of 
(4.05). They are recorded in the upper part of table 1, correct to 9 significant figures. After passing 
the last of the given values of r, namely 10, the "test function" pioe r / (p r Pr+i) was computed. This 
falls below the value 2 X 10~ 8 for the first time when r= 14. In accordance with the criterion of 
section 4 this is the value 4 to be assigned to N. The column of values y j r 4) was then generated by 
backward use of (4.04), beginning with 3^ 1 4 4) = 0. For r— 1(1)10 these are the wanted approximations 
to Er(l). 

To test the accuracy of the results, the computations were repeated for 7V=32 and 7V=34, 
using a time-sharing automatic computer and working to 36 floating binary figures, with an ex- 
ponent of 12 binary figures. As further checks the values for r= 1(1)5 were compared with the 
10-decimal^ values given in [9], and the values for r=10 and 11 computed from the expansions 
(6.01). The full results of these computations are not included here, but the digits in y r 14) which 
differ from those in the more accurate values of E r (l) are printed in italic type, and the difference 
4 14) between the two values is recorded, in units of the 9th decimal place, in the penultimate 
column of the upper part of table 1. As expected, this error does not exceed 2 X 10 -8 in absolute 
value within the wanted range r= 1(1)10. 



4 The more precise criterion of section 5 requires that the sum of />io<? r /(p r Pr+i ) and all subsequent values in this column be less than 2 X 10 8 . Inspection of 
table 1 indicates that this also is satisfied for N= 14. This would not be the case, however, if we reduced our error tolerance in y r from 2 X 10 -8 to 1 X 10~ 8 . 
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263-048 O - 67 - 5 



Table 1. Anger-Weber function E r (l) 



r 


Pr 


e r 


Pwe r 

PrPr+l 


yi4) 


IO9414) 


109£ 14 p r 




1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 




1 

2 
7 
40 
313 
3090 
36767 
5 11648 
81 49601 
1461 81170 
2.91547 380X109 
6.39942 424X101° 
1.53294 634X101 2 
3.97926 IO6XIO1 3 
1.11266 015 X10 15 
3.33400 119 X10 16 


-0.56865 6627 
0.70458 2918 
0.70458 2918 
9.61725 973 
9.61725 973 
4.08141 237 X10 2 
4.08141 237 x 10 2 
4.72213 396 X10 4 
4.72213 396 X10 4 
1.04236 156 X10 7 
1.04236 156 X10 7 
3.72252 015X10 9 
3.72252 015X10 9 
1.95553 042X1012 
1.95553 042X1012 
1.41863 843X101 5 


3.6X10- 3 
2.9X10-3 

5.5 X10- 6 
4.7 X 10- 6 
6.5X10-9 

5.6 x lO- 9 


0.43816 2436 
.17174 1955 
.24880 5382 
.04785 0795 
.13400 0978 
.01891 9443 
.09303 2343 
.01029 3811 
.07166 8637 
.00650 2177 
.05837 3706 
.00447 9865 
.04974 3054 
.00000 0000 


1 

12 

240 

5279 

12 6445 


1 

12 

240 

5279 

12 6444 




Pr 


e r 


e r 




r 


PrPr+i 




14 
15 

16 
17 
18 
19 


3.97926 106 X10 13 
1.11266 015 X10 15 
3.33400 119X101 6 
1.06576 772X101 8 
3.62027 625X10 19 
1.30223 368X1021 


1.95553 042 x lOi 2 
1.41863 843X1015 
1.41863 843 x 10 15 
1.35839 625X1018 
1.35839 625X1018 


4.41672 x IO-1 7 

3.82422 x IO-1 7 

0.00399 x IO-1 7 

.00352 x 10-17 

.00000 x 10-17 


£,4 = 8.24845x10-17 



It is of interest to apply the expansions of section 5 to this example. The necessary computa- 
tions for evaluating the expansion (5.02) are given in the lower part of table 1, and the values of 
l0 9 Ei 4 p r appear in the final column of the upper part of the same table. They agree with 10 9 e ( r 14) 
to within a unit. 



Example 2. Struve functions. 

The function H r (x) satisfies eq (2.01) with 



a r = c r = 1 , 



b r = 



2r 



dr- 



(hxY 



V^rfr-h 



For fixed x and large r we have ([8] , sec. 10.4) 



H P (*) 



V2 



x /ex\ r 

? rrr W 



(6.05) 



Since the complementary functions of the difference equation are again the Bessel functions 
J r (x) and Y r (x), the conditions of section 2 are satisfied, provided that Jo(x) ¥" 0. 

Let us evaluate H r (0.1) to 8 significant figures for all positive integer values of r such that 
|H r (0.1)| exceeds Va X 1(H°. 

The computations are shown in table 2. The value 



e = Ho(0. 1) =0.06359 12700 
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Table 2. Struve function H,((). 1 ) 

















10 9 € ( , 15) 




Pr 




dr 


?r 




y<15) 




r 


PrPr+l 


y r i5) 








0.63661 


9772 


0.06359 12700 








1 


1 


2.12206 


591 x LO" 2 


.04237 06109 




2.12065 160X LO" 3 





2 


20 


4.24413 


182 x LO" 4 


.03388 23473 




4.24211 125 x 10 •• 





3 


799 


6.06304 


546 x L0"« 


.02903 79740 




6.06080 029 x 10 < 


- 1 


4 


47920 


6.73671 


718 x l()-« 


.02580 97391 




6.73467 605 x 10" 





5 


38 32801 


6.12428 


835 x 10- 10 


.02346 24212 




6.12271 820 x 10" 


2 


6 


3832 32180 


4.71099 


104 X 10- 12 


.02165 70178 




4.70994 424 X LO" 13 


4 


7 


4.59840 288 x 10 10 


3.14066 


069 x 10- 14 


.02021 28155 




3.14004 492 x 10" 15 


4 


8 


6.43738 080 x LO 12 


1.84744 


746 x 10- 16 


.01902 35432 




1.84712 338 X 10~ 17 


- 1 


9 


1.02993 494 X LO 15 


9.72340 


768 X10- 19 


.01802 20955 




9.72186 442 x 10~ 20 





10 


1.85381 852 x LO" 


4.63019 


413 x 10- 21 


.01716 37415 




4.62952 313 x 10" 22 


4 


11 


3.70753 405 X ]()"> 


2.01312 


788 x 10-" 


.01641 73675 




2.01285 948 X 10~ 24 


4 


12 


8.1.5638 953 X 10 2 ' 


8.05251 


152 X10- 26 


.01576 05733 




8.05151 746 X 10" 27 


2 


13 


1.95749 641 X 10 24 


2.98241 


167 x 10- 28 


.01517 67673 


1.5 x 10-" 


2.98206 890 X 10" 29 


-1 


14 


5.08940 91 OX 10 26 


1.02841 


782 x 10- 30 


.01465 33634 


2.0 x 10- 58 


1.02829 540 X 10-» 


1 1520 


75 


1.42501 497 x 10*' 


3.31747 


684 x 10-" 


.01418 06180 


2.3 x 10-" 







16 


4.27499 402 X 10 3 ' 















was extracted from [9], and confirmed by evaluation of the expansion 



H r (*) = (**)'+» £ 



~r(. + gr(.+HJ) 



(6.06) 



The largest of the wanted values of r was determined by the criterion 

|e r /p r+1 |>|XlO-™ and |e r+1 /p r+2 |^ .} X 10-™; 

compare (5.03). This gave r= 13. Next, we have from (5.01), (5.02), and (5.03), 

^ A) = 2r£r±i g -v 
y r ' e r PnPn+i 

From table 2 we see that the right of this relation is an increasing function of r, hence /V is the 
least value for which 



e.\ 



p.\P\+\ 



(i xio- 8 ) 



gig 

PnPu 



From the entries in the column headed e r l(prPr+\) we see immediately that this gives yV= 15. 

The values of yV 5) , computed from (4.04), appear in the penultimate column of the table. For 
r ^ 13 they are the required approximations to H,(0.1). Again, more accurate values were obtained 
by automatic computation with a higher value of /V(26), and also by evaluation, of the expansion 
(6.06) for r= 1(1)15. In the final column the relative error e^V/,- 1 '^ is given in units of the 9th deci- 
mal place. As expected, it lies within the stipulated limit ¥2 X 10~ 8 in the required range. 

7. Propagation of Rounding Errors 

In addition to the truncation error e\: v) which has been analyzed in sections 4 and 5, the other 
possible sources of error in the final solution are the rounding errors introduced during the calcu- 
lations. Since the computing process is essentially the solution of a finite system of linear algebraic 
equations, the nature of the transmission of these errors is available from general theory [10] 
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chapter 4; [6], chapter 9. However, because of special features of the present problem, including 
the fact that in our form of elimination the absolute values of the multipliers are not bounded by 
unity, some comments on the effects of rounding errors may be helpful. 

Consider first the computation of the sequence p r . From the conditions p = 0, p\ = 1, we see 
that in terms of the fundamental solutions f r and g r of section 2 

Pr = (fogr - g(fr) I ifogl ~ g(fl ) , (7.01) 

the denominator here necessarily being nonzero since f r and g r are independent solutions of the 
difference equation. By hypothesis, /o ^ 0; therefore p r always contains a multiple of g r . And since 
frlgr~^0 as r— »o°, p r ultimately becomes proportional to g r when a fixed number of significant 
figures is maintained in the computations. 

Each rounding error in the formation of the p r can be regarded as introducing unwanted small 
multiples of f r and g r . Ultimately, the former dies out in comparison with the latter; the error is 
then propagated at the same rate as p r itself. Before this stage is attained, however, some loss of 
accuracy is possible. If the value of |/o| is unduly small compared with \gofilgi\, then from (7.01) 
we see that initially p r behaves like a multiple of f r . But the rounding errors are still propagated 
in proportion to g r , and this generally causes a steady loss of significant figures. The loss ceases 
when the term/ogr in (7.01) overtakes gqf r in magnitude, at which stage the computation becomes 
completely stable. 

It should be realized that this loss of accuracy is not attributable to the method of computation, 
but to the fact that, as a rule, the whole problem is ill-posed when \f \ is small compared with 
\gofrlgr\ for at least one value of r. For from (2.06) we see that 



3y r 


= 


frdk 
Mr 


> 


gM 
goYr 



(7.02) 

where 8y r is the change in y r consequent upon an arbitrary change dk in the value of k. Since 
|g>/gb| is generally large compared with \y r \ (see (2.04)), the relative error in y r is very sensitive 
to rounding errors in the given value of k. 

Examples 1 and 2 of section 6 would be ill-posed in this way if the chosen value of x were 
close to a zero of Jo(x), say x = 5.52. This would become apparent at the beginning of the compu- 
tations: the early p r would diminish in size, in contrast to the behavior they exhibit in tables 1 
and 2. 

The difficulty could be overcome in these and other examples by carrying out the computa- 
tion of k and p r to higher precision, and making the necessary prolongation of the recurrences until 
the criteria of sections 4 and 5 for terminating them are met. 

If the value of y\ can be found, however, a preferable alternative is to apply the algorithm of 
sections 3 and 4 with the given y\ as normalizing value, instead of yo — k. In effect, this means 
that the recurrences (4.05) are begun with pi = 0, p 2 =l, and e\ = y\. Subsequently the value of 
yo can be computed from y x and yi by a single backward application of (2.01). 

The other part of the elimination process is the computation of the right-hand sides e r . From 
(4.05) we see that in the inhomogeneous case instability could arise from this source if there were 
persistent heavy cancellation between a r e r -\ and d r p r - No naturally occurring examples of this 
phenomenon have been encountered so far, however. 

Lastly, we see from (4.04) that a rounding error introduced in y[ N) during the back-substitution 
is multiplied by the factor p r lp s when it is transmitted to y^ (r < s). Except when the problem is 
ill-posed, this factor decays with diminishing rata faster rate than y ( r jV) itself, because p r contains 
a substantial multiple of g r , and y r contains no multiple of this function. 

Summarizing this section, we have shown that unstable transmission of rounding errors can 
occur only when the original problem is ill-posed or when heavy cancellation takes place during 
the calculation of e r from the second of (4.05). 
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8. Comparison With the Algorithms of Miller and Shintani 

In section 4 we solved the set of eqs (3.01) and (3.02) by eliminating the variables in the order 
y/^, y/\ . . ., >^2 : we ma y ca ^ tms forward elimination. Suppose now that these variables 
are eliminated in the reverse order: backward elimination. The resulting set of pivotal equations 
can be expressed in the form 

tig^-idp^**** (r=N-2,N-3, . . .,<)), (8.01) 

where the quantities u\^ and i//" are defined by 

u$>, = l, u$> 2 = b.s-tla N -u ^l, = ds-Jax- u (8.02) 



an 



(I 



!#?,=* — — , ^A = — ~, r<JV-l . (8.03) 

(It should be observed that \Af* and ^: v) depend on /V as well as r, unlike the p r and e r of section 4.) 
The last of eqs (8.01) is used to begin the back-substitution. It yields 

where k is again the given value of y<>. Then y/*, y/*, . . ., y^l, may be computed by successive 
application of (8.01) with r = 1, 2, . . ., AT— 2. 

Thus the elimination process consists of constructing a sequence // ( ,: v) which satisfies the 
homogeneous form of the given difference equation (2.01) and the conditions a^ = 0, 0^=1. 
This is exactly the first stage of Miller's algorithm [1], [3]: the u { r N) are the so-called trial values. 
And in the homogeneous case, given by d r = 0, all the quantities i/ r ^ vanish, causing the formulas 
(8.01) for back-substitution to reduce to 

„(JV) „(N) „(V) ,.(N) L 



U y "> HVL', w 



4N) = _ L_ *LN) — — I IZI -J— J.N)= UN) 

Tr rf«, Tr ~ l d r N \ i££ 2 ' ' ' ity n Off* r ' 

This is the second stage of the Miller algorithm: k/u^ is the normalizing factor. 

Accordingly, in the homogeneous case the Miller recurrence algorithm can be regarded as the 
solution of the set of equations {3.01) and (3.02), with d r = 0, by backward elimination. In the inho- 
mogeneous case the solution by backward elimination, described above, can be regarded as a 
generalization of the Miller algorithm. 

Compared with the forward elimination process of section 4, the Miller algorithm suffers from 
the disadvantages that it does not determine automatically the correct value of TV, and if a second 
value of /V is used as a check on the adequacy of the original value, then the computations must 
begin afresh. The advantage of the Miller algorithm is that the process of back-substitution is less 
laborious; this advantage is restricted to the homogeneous case, however, and is offset if more than 
one trial value of /V has to be used. 

The method Shintani [5] 5 has developed for solving second-order linear difference equations 
in the homogeneous case consists of the use of the Miller algorithm preceded by two forward recur- 
rence processes to determine the optimum value of N. In the present notation, Shintani takes 
a, = 1 and d r = 0. His formulas for forward recurrence are given by ([5], Theorem 1) 



5 See also [4], section 4. 
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Pr+X (v) = b r+ iP r (v) ~ CrPr-X {v) , (8.04) 

where v=0 or 1, and 

P-!(0)=0, Po(0)=l; P (1) = 0, Pi(l) = L (8.05) 

It is easily verified, for example, that the quantities P r (0) appear when the forward elimination 
procedure is applied to eqs (4.01) with d r — and the multipliers chosen in such a way that the 
constant value — k is preserved on the right-hand sides. The resulting pivotal equations are in fact 

-Pr{Q)W> + c r Pr-i{OWr%= ~k. (r=l,2, . . ., TV - 1 ) . (8.06) 

In our notation 

Pr(0)=CiC 2 . . . C r pr + l- (8.07) 

From the computational standpoint, the evaluation of Shintani's sequence P r (0) may be com- 
pared with the evaluation of our sequence p r , the evaluation of his P r (l) with our e r , and the appli- 
cation of the Miller algorithm with our process of back-substitution. In the first stage the com- 
puting effort is identical, but in the second and third stages our method requires considerably 
less effort. 

9. More General Form of Normalizing Condition 

Let us consider now the solution of the difference eq (2.01) when (2.05) is replaced by the more 
general normalizing condition 

raoyoH-raiyi + 7tt 2 y2 + . • .=.£, (9.01) 

in which m , m u . . ., and k are given constants. We again suppose that the general solution of 
(2.01) has the form (2.02), but instead of the conditions imposed on/ r , g r , and h r in section 2, we 
assume that 

N 

2 rrirgr ^°° as 7V-» oo ? (9 02) 

and 

Jm r / r = F, ^rrtrh r = H, (9.03) 

r=0 r=0 

where F and H are finite, and F/0. Then (2.01) has a unique solution fulfilling (9.01). It is given by 

k-H r , , 
Yr=— j=r-fr + h r ; (9.04) 

compare (2.06). 

The obvious extension of the approach of section 3 is to solve the system of linear algebraic 
equations given by 

ory™! - b r y™ + c r y*>, = dr (r= 1, 2 N- 1), (9.05) 

N 

^mry^ N) =k, (9.06) 



and 



y/ )= 0. (9.07) 
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THEOREM 2. In addition to the other conditions of this section, assume that for all sufficiently 
large N the system of equations (9.05), (9.06), and (9.07) has a solution, that g N ¥* 0, and that 

^ 2 m* -+ 0, ^ 2 m r g r -> 0, (N -» oo). ( 9.08) 

Then if r £5 fixed and N — > 00, y ( r N) — * y r . 

This result may be established by expressing y^ in the form 

yW = A s fr + Bjvgr + fer. (9.09) 

Using (9.06) and (9.07), we find that 

iV / N \ /iV \7V 

Ajv 2 mrgr — gN ( JJ m r h r — k\ fs f 2 rn r h r — k\—h N J) "v/r 

#V 2 m rfr~f\ ^ "Wr giV 2 TUrfr-fN ^ m ^r 

r=0 r=0 r=() r=0 

In consequence of the assumed conditions, the denominators are asymptotic to Fgs as N — > oo. 
Hence A N -+ (k — H)/F. Next, the assumed conditions imply th&tf N lg N and h\/gx both tend to zero. 
Hence ##— ► 0. Comparison of (9.04) and (9.09) completes the proof. 

When the forward elimination process of section 4 is applied to eqs (9.05), (9.06), and (9.07), 
the following pivotal equations are obtained: 



Pr+itfn-prffli + qr ( 2 m ^s N) ) = e r (r=0, 1, . . ., N - 1 ) , (9.10) 

(compare (4.04)), where 

a\ai . . . a r , _ N 

fl>=l, q r = (r^l), (9.11) 

CiC 2 . . . C r 

Po = 0, pi = m , e = k, (9.12) 

and, if r ^ 1, 

b r p r — OrPr-l , grgr-j ~ d r p, (Ck , Q x 

Pr+i= — c h^ r m r , e r = . (9.13) 

Cr Cr 

In consequence of (9.07), the final equation of the form (9.10) reduces to 

P.v^ v :\=eiv-i. (9.14) 

This yields the value of y^l { ; thence 3^2 2 »^-s' • • •» ^o' V) ma y ^ e com P ute d from (9.10) by back- 
substitution. 

The value of /V may be determined in a similar way to that suggested in section 4. Suppose, 
for example, that all nonvanishing values of y r are needed to a fixed number of decimal places, 

123 



D, say — a common form of requirement with the present type of normalizing condition. Then TV 
is determined by the condition 



kv/p.v+i|< J xio-a, 
provided that \p r \ ^ \p s \ when r ^ N, and also 

\prlq r -\\ > \m r \, |m r +i|, . . . , |m. v |. 

Example 3. Bessel functions. 

Let us evaluate Jo (x) , J\ (x) , . . . , for x = 5 to 5 decimal places, by use of the relations 



Jr-l(x)-(2rlx) Jr(x)+J r+ l(x)=0, 



and 



Mx)+2Jt(x)+2J 4 (x) + . . . = 1. 

In the present notation, we have 

ar = c, = l, b r = 2r/x, d r = 0, 
mo=l, m\=mz = . . . = 0, m-i = m\ — . . . = 2, 
Accordingly, eqs (9.11) through (9.13) yield 

g r =l, e r =l, Po = 0, pi = l, 

and 



(9.15) 
(9.16) 

(9.17) 
(9.18) 



4=1. 



(9.19) 



Pr+i =b r pr — pr-i + rn r . 

Table 3 gives the values of p r correct to 6 significant figures. The criterion (9.15) suggests that 
TV be taken as the least value of r for which \p r +i | > 2 X 10 5 . This gives 6 N= 14. The column of 
values of y< r 14) is then generated upwards by use of (9.10), starting with y i l \ 4) = 0. These are the 
required approximations to A(5): their differences, e ( r 14) , from the true values are recorded in the 
final column in units of the 5th decimal place. The agreement is satisfactory. 

Table 3. Bessel function J r (5) 













13 




r 


b r 


m r 


Pr 


yiU) 


2 m ^ 14) 


1()6 € (14) 





0.0 


1 





-0.17758 




-2 


1 


0.4 





1 


-.32758 


1.17758 





2 


0.8 


2 


0.4 


.04655 


1.17758 


2 


3 


1.2 





1.32 


.36482 


1.08448 


1 


4 


1.6 


2 


1.184 


.39123 


1.08448 





5 


2.0 





2.5744 


.26114 


0.30202 





6 


2.4 


2 


3.9648 


.13105 


.30202 





7 


2.8 





8.94112 


.05338 


.03992 





8 


3.2 


2 


21.0703 


.01841 


.03992 





9 


3.6 





60.4838 


.00552 


.00310 





10 


4.0 


2 


196.671 


.00147 


.00310 





11 


4.4 





728.200 


.00035 


.00016 





12 


4.8 


2 


2007.41 


.00008 


.00016 





13 


5.2 





13709.4 


.00001 


.00000 


1 


14 


5.6 


2 


68281.5 


.00000 


.00000 





15 






368669. 









fi When the condition (9. 16) is violated — as it is in this example near the beginning of the range — it would be safer in practice to take N slightly higher than the 
value predicted by (9.15). 
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It may be noted that estimates of the optimum value of TV for generating Bessel functions from i 
(9.17) and (9.18) by Miller's algorithm have been computed by Makinouchi [11] for x = 0.01 (.01) 
0.1(. 1)1(1)10(10)100 and precisions of 9, 10, 18, 20, and 30 significant figures. These values were 
obtained by use of the asymptotic approximations (6.04) above. In constructing a program for 
generating the J r (x) for arbitrary x and arbitrary precision, however, it would be simpler to deter- 
mine the optimum N by use of (9.15) (or (4.11)). The resulting gain would tend to offset the extra 
effort needed in applying the back-substitution relation (9.10) compared with the normalizing of 
the trial values in the Miller algorithm. 



10. Bounds for the Truncation Error 

In order to obtain strict bounds for the truncation error associated with the algorithm of section 
9, we proceed as in section 5. Write, temporarily, 

7) r = ^-^\ 



Then from (9.10) we obtain 

Therefore 

tor|*Pr(for+i| + for+l| + - • -+M) fr < N) * 

where p, is the greater of 

p r — q r m r +i 

id | - i — 



Pr+l 



lPr+1 



sup \m r +s\. 

2s£.V=£oc 



(10.01) 



(10.02) 



(10.03) 



(10.04) 



Equations (9.07), (9.14), and (10.01) yield T7# = ejv/pAM-i. From this result and (10.03) we may 
verify that 

\r) N -i\ ^pv-ikv/p,v+i|< (10.05) 

and thence by induction that 

|l) r | ^pr(l+Pr+l)(l+Pr+ 2 ) ■ ■ • (1 + ptf-l) MpAH-l| (r ^ N ~ 2) . (10.06) 

The left-hand side o/ the last relation is \y^ +1) — yM|. Replacing N by N+ 1, /V+2, . . . , in turn 
and summing, and applying Theorem 2, we find that 



I^I^PrU + Pr+OU+Pr^) . . . {l+fto-i)E N (r^/V~2), 



where 



and 



6 0v) = yr _ r (Ar) ? 



E N 



es 



P\+\ 



U+p.v) 



ey+i 



PN+2 



+ (l + p N )(l + p N+t ) 



eiv+s 



PN+3 



+ . . ., 



provided that the last series converges. Similarly 



(10.07) 
(10.08) 

(10.09) 
(10.10) 
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The results (10.07) and (10.10) are strict bounds for the truncation error, in contrast to the 
expansion of section 5 which is exact (for the algorithm of sections 3 and 4). Often the bound (10.07) 
is a considerable overestimate. 7 Thus in Example 3, the right-hand side of (10.07) or (10.10) has the 
following values for N = 14, in units of the 5th decimal place: 

568, 237, 29, 12, 3, 1, 1, then zero for r= 7, 8, . . . , 14. 

In consequence, if N is determined by the criterion that for each required value of r the right-hand 
sides of (10.07) and (10.10) must not exceed the specified tolerance in y r , then the resulting value 
is perfectly safe but often unnecessarily high. Applied to Example 3, this criterion yields A^=18, 
compared with the value 14 which we used and found to be quite adequate. 

In the next section we give an alternative formulation of the algorithm of section 9. Although 
perhaps less elegant, it generally yields a sharper assessment of the truncation error than that 
of this section. 

1 1 . Alternative Method for the General Normalizing Condition 

The algorithm of sections 3 and 4 can be applied to the problem of section 9 in the following 
way. First, we construct a solution f r of the homogeneous form of the given equation (2.01). The 
choice of this solution is arbitrary, provided that the first of (2.03) is satisfied. Then by means of 
an additional back-substitution we construct an arbitrary solution h r of (2.01) itself. The required 
solution y r may then be computed from (9.04), in which k is defined by (9.01), and F, H by (9.03). 
In the case when the given difference eq (2.01) is itself homogeneous, only the solution f r need be 
computed, and (9.04) reduces to 

y r = {klF)f r . (11.01) 

The simplest choice of the normalizing conditions needed for constructing/ r and h r is given by 

/o=l, h = 0. (11.02) 

The first of these may be an inconvenient or even impossible condition, however; in this event we 
may follow the suggestion given in section 7 and use instead 

/i = l, h, = 0. (11.03) 

To assess the truncation error in the final solution y r , let <p^ and 6[ N) be the truncation errors 
in the approximations / ( r v) and h\. N) to f r and h r ; thus 

fr=f<r N) + <P ( r N \ h r = h^+6["K (11.04) 

Bounds for <£ ( ,. v) and 0\: x) are computable from the expansions of section 5. From (9.04) we have 

yr = k ~ H *~ T * (fr N) + <P ( r N) ) + W/° + W (r= 0, 1, . . . , TV), (l 1.05) 

r a -r <T\ 

where F\, H\ are the computed quantities 

F .V= 2' m rfr N K H,= 2 mrhM, (1 1.06) 



'This can be traced to the fact that over most of the range the second of the two quantities (10.04) is usually very much smaller than the first. 
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and o\v, T\ assessable errors 5 



°".v = E W ^ ( / ;V) + 2 ,7lr ^ r ' (11.07) 



T.v=2 ^r»V V) + £ "*r/l,- (11.08) 



If F\ ¥^ 0, then to the first order of small quantities the truncation error in the formula 

-r 
F N 

is composed of three parts: 



rr-^k/vn+Ay" < 11091 



i^n, _{ T , +2 *!£*J)£:. „»,, (rSN) , (1U0) 

In the homogeneous case they reduce to two: 

i-<fV>, -TTJV\ (r«AT). (11.11) 

The first of (11.11) is the normalized multiple of the truncation error in the formula f r ^f { r N) ; the 
second of (11.11) is a fixed relative error arising from the approximate representation of the nor- 
malizing factor k/F by ///\ v . 

The equivalence of the method of this section to the algorithm of section 9 can be seen from 
the fact that the function on the right of (11.09) is exactly the solution of the set of eqs (9.05), 
(9.06), and (9.07). 

Example 4. J) 

Let us compute to 5 decimal places the solution of the homogeneous equation 

(2r-l)y r _,-12ry r +(2r+l)y r+1 = 0, (11.12) 

satisfying the condition 

Jyo + yi + y 2 + y 3 +. • ■ = !. (11.13) 

In the notation of sections 2 and 9 we have 

a r = 2r-l, 6 r =<12r, c r =2r+l, d r = 0, m = %, m r =\ (r>0), k=\. 

The computations are shown in table 4. Values of p r were generated from p — 0, pi = l, and (11.12) 
when r> 1, correct to 6 significant figures. With eo= 1 (compare the first of (11.02)), we find from 
the second of (4.05) that e r = l/(2r+ 1). The least value of r for which e r /p r +i < % X 10~ 5 is 7; in 



8 See also [7J. 

9 [3], section 5. 

10 An alternative way of estimating N in this particular example is to observe that, except for small r, the wanted solution of (11.12) behaves roughly like Ak- r , 
where A is a constant and X* — 6X + 1 = 0, giving X*^ 3 + V»=5.8. Assuming that A is of order unity, as is reasonable in view of the condition (11.13), we find that 
7V^5/logio X%7. This method of estimation is somewhat less certain than the one we have used, and it is not universally applicable. 
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accordance with (4.12) this is the value 10 to ascribe to N. The back-substitution process for the 
determination of/J. 7) is given byf\ 7) = 0, and 



Pr+l/™ = Pr/™+e r 



(r = 6, 5, . . .,0); 



compare (4.04). Division of / r 7) by F 7 = 0.599069, computed from the first of (11.06), yields the 
wanted approximations y r 7) to y r . 



Table 4 



r 


Pr 


e r 


JV ] 


y<) 








1.000000 


1.000000 


1.66926 


1 


1 


0.333333 


0.086107 


0.14373 


2 


4 


.200000 


.011094 


.01852 


3 


18.6 


.142857 


.001587 


.00265 


4 


92.8 


.111111 


.000238 


.00040 


5 


480.467 


.090909 


.000037 


.00006 


6 


2544.80 


.076923 


.000006 


.00001 


7 


13687.7 


.066667 


.000000 


.00000 


8 


74445.6 









The example is now complete, but it is of interest to illustrate the error analysis of this section. 
Accordingly, the whole calculation was repeated twice, keeping four extra significant figures 
throughout. In the first repetition the same value N—l was used. In the second repetition a new 
N was determined by the condition \eNlpN+\\ < § X 10~ 9 ; this gave A^= 12. 

The results appear in table 5. The column headed 10 9 €{. 7) gives the difference of Wy ( r 7) from 
the more accurate values 10 9 yp 2) . The next columns give 10Vr 7) /^7 and — l0 9 <r 7 f^lF*; the value 
of <£> ( r 7) was obtained by subtracting / ( r 7) from /V 2) , and cr 7 computed from (11.07), using the values 
of /V 2) for f r when r ^ 8. As expected, the values of 10 9 e ( r 7) are in good agreement with the sum of 
the entries on the same row in the following two columns. 



Table 5 



r 


fin) 


y r i2) 


fi 7) 


y?) 


io 9 4 7) 


6P 

F 7 


-10»%f 
Fi 





1.00000 00000 


1.66925 3684 


1.00000 00000 


1.66925 7339 


-3655 





-3655 


1 


0.08610 68379 


0.14373 4156 


0.08610 68378 


0.14373 4471 


-315 





-315 


2 


.01109 40183 


.01851 8731 


.01109 40180 


.01851 8771 


-40 


1 


-41 


3 


.00158 71852 


.00264 9415 


.00158 71839 


.00264 9418 


-3 


2 


-6 


4 


.00023 83677 


.00039 7896 


.00023 83614 


.00039 7887 


9 


11 


-1 


5 


.00003 68169 


.00006 1457 


.00003 67845 


.00006 1403 


54 


54 





6 


.00000 57914 


.00000 9667 


.00000 56199 


.00000 9381 


286 


286 





7 


.00000 09228 


.00000 1540 


.00000 00000 


.00000 0000 


1540 


1540 





8 


.00000 01485 
-00000 00241 


.00000 0248 
.00000 0040 












9 










10 


.00000 00039 


.00000 0007 


F 


'7 = 0.59906 88055 


o- 7 = 0.00000 13118 




11 


.00000 00006 


.00000 0001 


F 


2 = 0.59907 01173 






12 


.00000 00000 


.00000 0000 











12. Summary 

In this paper we have described a new algorithm for computing the solution y r of any second- 
order linear difference equation, homogeneous or inhomogeneous, which is applicable when simple 
forward recurrence (and possibly also backward recurrence) cannot be used because of instability. 
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In the first part (sees. 2-8) we considered the case in which the wanted solution y r has a 
specified value at the beginning of the range r=0, and an appropriate convergence condition as 
r ^oc. In this case the algorithm is based on the solution of a finite number, N 9 of simultaneous 
linear algebraic equations of tridiagonal form by forward elimination. As N^>°° the solution 
>i Ar) of these equations converges to y r (sec. 3). In sections 4 and 5 it was shown that during the 
process of computing y { r N) the minimum value of A^ necessary to achieve specified tolerance in 
\jr~ yV^I emerges automatically. Analyses of the truncation error and of the propagation of 
rounding errors were made in sections 5 and 7. The former leads to a convergent series expansion 
for y r ; the latter shows that the method of computation is quite stable, unless the problem itself 
is ill-posed. Numerical examples (sec. 6) illustrated the algorithm and confirmed the error analyses. 

In section 8 it was shown that the well-known algorithm of J. C. P. Miller for the homogeneous 
case can be regarded as the computation of y^ r N) by backward elimination, taking a guessed value 
of N. It was also shown that the recent extension of Miller's algorithm by Shintani is related to the 
process of forward elimination. 

In the second part of the paper (sees. 9-11) a more general form of normalizing condition for 
y r was considered. An extended form of the algorithm was developed in section 9 and applied to a 
numerical example in the same section. In section 10 bounds for the truncation error were given 
and discussed. In the concluding section (sec. 11) it was shown that the more general problem can 
also be solved by application of the original algorithm of sections 3 and 4. 

It is hoped that the results of this paper will prove to be of considerable usefulness in the 
computation of special functions from recurrence relations, in the solution of ordinary differential 
equations in Chebyshev series by Clenshaw's method, and in the solution of the discretized form 
of boundary-value problems in ordinary differential equations when one boundary is at infinity. 
In the last two connections, it may be possible to extend the present approach to difference equa- 
tions of order higher than the second. 
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